#Assessing potential influence of unusual values in final data set

hist(final_data$mg.C.per.g.soil_FLF) #very high concentration of FLF C in one sample

hist(final_data$mg.C.per.g.soil_OLF)

hist(final_data$mg.C.per.g.soil_HF)

hist(final_data$mg.N.per.g.soil_FLF) #very high concentration of FLF N in one sample

hist(final_data$mg.N.per.g.soil_OLF)

hist(final_data$mg.N.per.g.soil_HF)

hist(final_data$proportion.C_FLF) 

hist(final_data$proportion.C_OLF)

hist(final_data$proportion.C_HF)

hist(final_data$proportion.N_FLF) 

hist(final_data$proportion.N_OLF)

hist(final_data$proportion.N_HF)

#removing HARV_037 from final dataset


mean_FLF_conc.C=mean(final_data$mg.C.per.g.soil_FLF) #2.09
std_dev_FLF_conc.C=sd(final_data$mg.C.per.g.soil_FLF) #1.98
#HARV_037 is 47.33, greater than 20 standard deviations above mean

mean_FLF_conc.N=mean(final_data$mg.N.per.g.soil_FLF) #0.106
std_dev_FLF_conc.N=sd(final_data$mg.N.per.g.soil_FLF) #0.205
#HARV_037 has 1.42 mg FLF N per g soil, greater than 6 standard deviations above mean

final_data=final_data %>% 
  filter(plotID != "HARV_037") 

Site Characteristics

Site MAT (C) MAP(mm) Dominant tree species Soil % clay (Mean +/- SD) Soil Fe (ox-extractable)
LENO 18.1 1386 American Sweetgum, American Hornbeam, Possumhaw 44.7 +- 13.8 1.0 +- 0.4
DELA 17.6 1372 Water oak, Red maple, Sugarberry 37.0 +- 6.3 0.88 +- 0.2
ORNL 14.4 1340 Red maple, Sour gum, Chestnut oak 22.0 +- 12.7 0.13 +- 0.05
SERC 13.6 1075 Tulip poplar, American Beech, American Sweetgum 17.4 +- 6.2 0.27 +- 0.1
HARV 7.4 1199 Eastern Hemlock, Northern Red Oak, American Beech 5.5 +- 1.9 0.38 +- 0.1
BART 6.2 1325 American Beech, Eastern Hemlock, Red maple 3.6 +- 0.7 0.76 +- 0.5
TREE 4.8 797 Sugar maple, Red maple, Gray alder 5.4 +- 1.1 0.34 +- 0.09
  proportion MAOM C [MAOM C] proportion MAOM N [MAOM N]
Predictors Estimates std. Error std. Beta standardized std. Error p df Estimates std. Error std. Beta standardized std. Error p df Estimates std. Error std. Beta standardized std. Error p std. p df Estimates std. Error std. Beta standardized std. Error p df
Intercept 0.80 0.04 0.00 0.15 <0.001 41.00 30.15 7.53 -0.00 0.24 <0.001 41.00 0.88 0.04 0.02 0.11 <0.001 0.847 40.00 1.18 0.34 -0.00 0.18 0.001 41.00
Mycorrhizal dominance -0.10 0.03 -0.42 0.13 0.002 41.00 -1.96 3.05 -0.07 0.10 0.524 41.00 -0.19 0.07 -0.23 0.12 0.010 0.061 40.00 -0.14 0.19 -0.08 0.11 0.453 41.00
CDI 0.15 1.23 0.02 0.18 0.905 41.00 -583.93 212.12 -0.69 0.25 0.009 41.00 0.72 1.16 0.50 0.13 0.536 0.001 40.00 -15.46 9.69 -0.32 0.20 0.118 41.00
FeOx 0.01 0.00 0.50 0.16 0.003 41.00 1.36 0.33 0.62 0.15 <0.001 41.00 0.00 0.00 0.09 0.14 0.517 0.517 40.00 0.11 0.02 0.88 0.15 <0.001 41.00
Mycorrhizal dominance *
CDI
4.33 1.99 0.24 0.11 0.036 0.036 40.00
Random Effects
σ2 0.00 36.91 0.00 0.14
τ00 0.00 Site 27.96 Site 0.00 Site 0.04 Site
ICC 0.09 0.43   0.24
N 7 Site 7 Site 7 Site 7 Site
Observations 47 47 47 47
Marginal R2 / Conditional R2 0.330 / 0.391 0.359 / 0.635 0.429 / NA 0.474 / 0.602
  proportion oPOM C [oPOM C] proportion oPOM N [oPOM N]
Predictors Estimates std. Error std. Beta standardized std. Error p df Estimates std. Error std. Beta standardized std. Error p df Estimates std. Error std. Beta standardized std. Error p std. p df Estimates std. Error std. Beta standardized std. Error p df
Intercept 0.06 0.03 0.00 0.18 0.039 41.00 3.55 1.33 0.00 0.19 0.011 41.00 0.03 0.03 -0.03 0.13 0.227 0.847 40.00 0.13 0.04 0.00 0.20 0.005 41.00
Mycorrhizal dominance 0.05 0.02 0.33 0.15 0.027 41.00 1.25 0.80 0.20 0.13 0.127 41.00 0.11 0.04 0.18 0.14 0.016 0.191 40.00 0.02 0.03 0.09 0.13 0.494 41.00
CDI 0.48 0.84 0.12 0.20 0.574 41.00 -79.58 37.57 -0.44 0.21 0.040 41.00 0.10 0.69 -0.40 0.16 0.881 0.014 40.00 -2.62 1.20 -0.48 0.22 0.035 41.00
FeOx -0.00 0.00 -0.17 0.18 0.347 41.00 0.12 0.08 0.25 0.17 0.148 41.00 0.00 0.00 0.09 0.16 0.590 0.590 40.00 0.00 0.00 0.30 0.18 0.099 41.00
Mycorrhizal dominance *
CDI
-2.56 1.18 -0.28 0.13 0.036 0.036 40.00
Random Effects
σ2 0.00 2.63 0.00 0.00
τ00 0.00 Site 0.56 Site 0.00 Site 0.00 Site
ICC 0.10 0.18   0.19
N 7 Site 7 Site 7 Site 7 Site
Observations 47 47 47 47
Marginal R2 / Conditional R2 0.101 / 0.190 0.215 / 0.353 0.246 / NA 0.189 / 0.346
  proportion.C_FLF mg.C.per.g.soil_FLF proportion.N_FLF mg.N.per.g.soil_FLF
Predictors Estimates std. Error std. Beta standardized std. Error p df Estimates std. Error std. Beta standardized std. Error p df Estimates std. Error std. Beta standardized std. Error p df Estimates std. Error std. Beta standardized std. Error p df
(Intercept) 0.137 0.029 -0.001 0.132 <0.001 41.000 4.720 1.502 -0.009 0.216 0.003 41.000 0.119 0.021 0.000 0.122 <0.001 41.000 0.174 0.051 -0.009 0.211 0.001 41.000
E 0.051 0.023 0.280 0.127 0.033 41.000 1.011 0.745 0.162 0.119 0.182 41.000 0.027 0.018 0.198 0.128 0.129 41.000 0.019 0.027 0.091 0.124 0.471 41.000
CDI -0.715 0.809 -0.136 0.154 0.382 41.000 -87.344 42.376 -0.480 0.233 0.046 41.000 -1.610 0.577 -0.400 0.143 0.008 41.000 -3.068 1.429 -0.491 0.229 0.038 41.000
FeOx -0.007 0.002 -0.498 0.150 0.002 41.000 -0.021 0.078 -0.044 0.164 0.791 41.000 -0.003 0.002 -0.241 0.145 0.105 41.000 0.000 0.003 0.008 0.170 0.965 41.000
Random Effects
σ2 0.002 2.223 0.001 0.003
τ00 0.000 Site 0.951 Site 0.000 Site 0.001 Site
ICC 0.031 0.300   0.261
N 7 Site 7 Site 7 Site 7 Site
Observations 47 47 47 47
Marginal R2 / Conditional R2 0.346 / 0.366 0.265 / 0.485 0.333 / NA 0.233 / 0.433
#what % of total C and N was MAOM C in each site?
pct.MAOM.C=final_data %>% 
  group_by(Site) %>% 
  mutate(mean.C.pro.maom=mean(proportion.C_HF),
         se.pro.maom=sd(proportion.C_HF)/sqrt(46)) %>% 
  dplyr::select(plotID,mean.C.pro.maom, se.pro.maom)
## Adding missing grouping variables: `Site`
pct.MAOM.N=final_data %>% 
    group_by(Site) %>% 
  mutate(mean.N.pro.maom=mean(proportion.N_HF),
         se.pro.maom.N=sd(proportion.N_HF)/sqrt(46)) %>% 
   dplyr::select(plotID,mean.N.pro.maom, se.pro.maom.N)
## Adding missing grouping variables: `Site`
# what % of total C and N was MAOM at the lowest %ECM and at the highest %ECM in each site?

end_members_maom_prop=final_data %>% 
  select(Site,Plot.x,E, proportion.C_HF, proportion.N_HF) %>%
  arrange(Site,E) %>% 
   group_by(Site) %>% 
  mutate(myc_rank=order(E, decreasing=T)) %>% 
  mutate(max_A=max(myc_rank)) %>% 
  filter(myc_rank==max_A) %>% 
  ungroup() %>% 
  mutate(mean_prop_C=mean(proportion.C_HF),
         mean_prop_N=mean(proportion.N_HF))

 #At max ECM dominance, avg. proportion C in MAOM was 0.75831, N was 0.85565

#At max AM dominance, avg proportion C in MAOM was 0.83787, N was 0.89959

#from AM to ECM, MAOM C prop dropped 7.956%, N dropped 4.39%

concC.fe.hf=ggplot(final_data, aes(x=FeOx, y=mg.C.per.g.soil_HF, color=Site))+
   geom_smooth(method="lm",size=2, se=F)+
  geom_point()+
   scale_colour_viridis_d(option="C")+
  labs(x=" FeOx (mg/g soil)", y="[MAOM C] (mg/g soil)", colour="Site")+
  theme_cowplot()+
  theme(legend.title=element_text(hjust=0.5, size=16), legend.text=element_text(size=15), axis.text=element_text(size=15), axis.title=element_text(size=19))

propC.fe.hf=ggplot(final_data, aes(x=FeOx, y=proportion.C_HF, color=Site))+
   geom_smooth(method="lm",size=2, se=F)+
  geom_point()+
   scale_colour_viridis_d(option="C")+
  labs(x=" FeOx (mg/g soil)", y=expression("proportion C"[MAOM]), colour="Site")+
  theme_cowplot()+
  theme(legend.title=element_text(hjust=0.5, size=16), legend.text=element_text(size=15), axis.text=element_text(size=15), axis.title=element_text(size=19))


concN.fe.hf=ggplot(final_data, aes(x=FeOx, y=mg.N.per.g.soil_HF, color=Site))+
   geom_smooth(method="lm",size=2, se=F)+
  geom_point()+
   scale_colour_viridis_d(option="C")+
  labs(x=" FeOx (mg/g soil)", y="[MAOM N] (mg/g soil)", colour="Site")+
  theme_cowplot()+
  theme(legend.title=element_text(hjust=0.5, size=16), legend.text=element_text(size=15), axis.text=element_text(size=15), axis.title=element_text(size=19))

ggarrange(concC.fe.hf,concN.fe.hf, nrow=1, ncol=2, common.legend=T, legend="right")
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

ggsave("hf.conc.CN.fe.pdf",path="figures/", width= 20, height= 10, units="cm")

#what are the slopes of these lines?
fe_conchf_slopes_C=coef(lmList(mg.C.per.g.soil_HF~FeOx|Site , data = final_data))[2]
fe.mod1.slopes= as.data.frame(forests)%>% 
  arrange(forests) %>% 
  cbind(fe_conchf_slopes_C$FeOx) %>% 
  summarize(mean_slope=mean(fe_conchf_slopes_C$FeOx))

fe_conchf_slopes_N=coef(lmList(mg.N.per.g.soil_HF~FeOx|Site , data = final_data))[2]
fe.mod2.slopes= as.data.frame(forests)%>% 
  arrange(forests) %>% 
  cbind(fe_conchf_slopes_N$FeOx) %>% 
  summarize(mean_slope=mean(fe_conchf_slopes_N$FeOx))
##         1         2 
## -26.40569 -26.45323
  d13C d14C
Predictors Estimates std. Error std. Beta standardized std. Error p Estimates std. Error std. Beta standardized std. Error p
(Intercept) -26.25 0.24 -0.00 0.14 <0.001 38.71 51.93 -0.01 0.40 0.461
FeOx 0.07 0.02 0.64 0.20 0.002 1.89 1.39 0.19 0.14 0.185
CDI -17.29 7.45 -0.45 0.19 0.026 -1579.98 1489.10 -0.44 0.42 0.296
E -0.03 0.21 -0.02 0.15 0.904 -12.25 13.00 -0.10 0.10 0.353
Random Effects
σ2 0.16 564.42
τ00 0.00 Site 1567.09 Site
ICC   0.74
N 6 Site 6 Site
Observations 40 40
Marginal R2 / Conditional R2 0.230 / NA 0.089 / 0.759
  d14C d14C d13C d13C
Predictors Estimates std. Error std. Beta standardized std. Error p Estimates std. Error std. Beta standardized std. Error p Estimates std. Error std. Beta standardized std. Error p Estimates std. Error std. Beta standardized std. Error p
(Intercept) -12.86 14.42 0.00 0.16 0.378 -56.20 74.97 0.00 0.16 0.458 -26.42 0.16 -0.00 0.16 <0.001 -27.08 0.80 -0.00 0.16 <0.001
mg C per g soil HF 0.21 0.74 0.05 0.16 0.775 -0.00 0.01 -0.01 0.16 0.973
proportion C HF 58.01 92.10 0.10 0.16 0.533 0.81 0.99 0.13 0.16 0.419
Observations 40 40 40 40
R2 / R2 adjusted 0.002 / -0.024 0.010 / -0.016 0.000 / -0.026 0.017 / -0.009
  d14C d14C d14C d14C d14C d14C
Predictors Estimates std. Error std. Beta standardized std. Error p Estimates std. Error std. Beta standardized std. Error p Estimates std. Error std. Beta standardized std. Error p Estimates std. Error std. Beta standardized std. Error p Estimates std. Error std. Beta standardized std. Error p Estimates std. Error std. Beta standardized std. Error p
(Intercept) -96.48 51.54 -0.00 0.42 0.135 -27.34 22.13 0.00 0.43 0.284 43.34 26.87 -0.00 0.33 0.182 48.99 65.67 -0.00 0.46 0.510 -118.57 64.15 0.00 0.35 0.138 -7.89 32.31 -0.00 0.43 0.823
FeOx 4.69 5.25 0.41 0.45 0.422 0.52 2.99 0.13 0.74 0.870 5.40 6.35 0.43 0.50 0.443 -10.42 11.74 -0.48 0.54 0.440 6.62 6.00 0.66 0.59 0.332 7.68 9.80 0.44 0.56 0.490
E -22.78 49.88 -0.21 0.45 0.672 17.43 45.33 0.28 0.74 0.720 -17.01 26.68 -0.32 0.50 0.558 -26.97 46.72 -0.31 0.54 0.604 115.10 68.60 1.00 0.59 0.169 -36.78 29.43 -0.70 0.56 0.300
Observations 7 7 7 6 7 6
R2 / R2 adjusted 0.187 / -0.219 0.156 / -0.267 0.480 / 0.219 0.226 / -0.290 0.420 / 0.130 0.344 / -0.093
  d13C d13C d13C d13C d13C d13C
Predictors Estimates std. Error std. Beta standardized std. Error p Estimates std. Error std. Beta standardized std. Error p Estimates std. Error std. Beta standardized std. Error p Estimates std. Error std. Beta standardized std. Error p Estimates std. Error std. Beta standardized std. Error p Estimates std. Error std. Beta standardized std. Error p
(Intercept) -26.31 0.40 -0.00 0.44 <0.001 -27.23 0.46 0.00 0.32 <0.001 -27.17 0.46 -0.00 0.16 <0.001 -26.87 0.67 0.00 0.51 <0.001 -26.48 0.50 0.00 0.29 <0.001 -27.25 0.50 -0.00 0.31 <0.001
FeOx 0.00 0.04 0.01 0.48 0.981 0.05 0.06 0.42 0.56 0.489 0.30 0.11 0.67 0.25 0.053 0.05 0.12 0.24 0.59 0.713 -0.03 0.05 -0.30 0.49 0.571 0.28 0.15 0.73 0.40 0.167
E -0.23 0.38 -0.29 0.48 0.579 0.57 0.94 0.34 0.56 0.577 -0.63 0.46 -0.34 0.25 0.239 -0.01 0.48 -0.02 0.59 0.978 0.57 0.53 0.53 0.49 0.341 0.15 0.46 0.14 0.40 0.759
Observations 7 7 7 6 7 6
R2 / R2 adjusted 0.084 / -0.374 0.521 / 0.281 0.877 / 0.815 0.061 / -0.565 0.611 / 0.416 0.664 / 0.440

ecm_order=received_samples_plot %>% 
  select(plotID, order) %>% 
  distinct()
## Adding missing grouping variables: `Site`
pctAngio=ggplot(plotPctAngio %>%
                  rename("Angiosperm"="fractionAngiosperm") %>% 
                  mutate(Gymnosperm=1-Angiosperm,
                         plotID=as.character(plotID))%>% 
                           pivot_longer(c(Angiosperm, Gymnosperm), names_to="gymAng", values_to="dominance") %>% 
                  left_join(ecm_order) , aes(x=order, y=dominance, fill=gymAng))+
  geom_bar(position="stack", stat="identity", colour="black")+
  scale_fill_manual( values=c("gray90", "gray20"))+
  xlab("Study Plot")+
  ylab("Proportion of Basal Area")+
  facet_wrap(~Site, scales="free_x")+
  theme_cowplot()+
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank(),strip.background = element_rect( fill="white"),legend.title=element_blank())
## Joining, by = "plotID"
ggsave("pctAngioGymno.jpg", path="figures/", width= 20, height= 16, units="cm")

BAbyMYC=ggplot(final_data %>% 
                 left_join(ecm_order), aes(x=order, y=plotBA))+
  geom_bar( stat="identity", colour="black")+
  xlab("Study plots in order of increasing ECM dominance")+
  ylab("Basal Area")+
  facet_wrap(~Site, scales="free_x")+
  theme_cowplot()+
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank(),strip.background = element_rect( fill="white"),legend.title=element_blank())
## Joining, by = c("plotID", "Site")
ggsave("BAmyc.jpg", path="figures/", width= 20, height= 16, units="cm")


ggplot(final_data, aes(x=fractionAngiosperm, y=proportion.C_HF, colour=Site))+
   geom_smooth(method="lm",size=2, se=F)+
   geom_point(size=3, alpha=0.5)+
  labs(x=" Angiosperm dominance ", y="MAOM C (proportion)", colour="Site")+
  theme_cowplot()+
  theme(legend.title=element_text(hjust=0.5, size=16), legend.text=element_text(size=15), axis.text=element_text(size=15), axis.title=element_text(size=19))
## `geom_smooth()` using formula 'y ~ x'

ggsave("maom.ang.prop.jpg", path="figures/", width= 19, height= 16, units="cm")
## `geom_smooth()` using formula 'y ~ x'
ggplot(final_data, aes(x=fractionAngiosperm, y=mg.C.per.g.soil_HF, colour=Site))+
   geom_smooth(method="lm",size=2, se=F)+
   geom_point(size=3, alpha=0.5)+
  labs(x=" Angiosperm dominance ", y="[MAOM C] (mg/g soil)", colour="Site")+
  theme_cowplot()+
  theme(legend.title=element_text(hjust=0.5, size=16), legend.text=element_text(size=15), axis.text=element_text(size=15), axis.title=element_text(size=19))
## `geom_smooth()` using formula 'y ~ x'

ggsave("maom.ang.conc.jpg", path="figures/", width= 19, height= 16, units="cm")
## `geom_smooth()` using formula 'y ~ x'
  mg.C.per.g.soil_HF mg.N.per.g.soil_HF
Predictors Estimates std. Error std. Beta standardized std. Error p Estimates std. Error std. Beta standardized std. Error p
(Intercept) 15.364 3.487 0.017 0.285 <0.001 1.245 0.185 -0.002 0.204 <0.001
plotBA 0.539 1.329 0.055 0.135 0.687 -0.042 0.087 -0.074 0.154 0.632
Random Effects
σ2 51.014 0.251
τ00 40.908 Site 0.044 Site
ICC 0.445 0.150
N 7 Site 7 Site
Observations 47 47
Marginal R2 / Conditional R2 0.003 / 0.447 0.005 / 0.155
#Relationships between site-level climate decomposition index (CDI) and a) soil oxalate-extractable iron content, and b) total tree basal area in our study plots.

basal_area_cdi=ggplot(final_data, aes(x=CDI, y=plotBA))+
   geom_point()+
  stat_summary(fun= mean, fun.min=mean, fun.max=mean, geom="crossbar", width=0.0009, position="dodge")+
  stat_summary(fun.data = mean_se, geom = "errorbar", width=0.0005, position = position_dodge(width = 0.0009))+
                theme_cowplot()+
  xlab("Climate Decomposition Index")+
  ylab("Total Basal Area"
  )


feox_cdi=ggplot(final_data, aes(x=CDI, y=FeOx))+
   geom_point()+
  stat_summary(fun= mean, fun.min=mean, fun.max=mean, geom="crossbar", width=0.0009, position="dodge")+
  stat_summary(fun.data = mean_se, geom = "errorbar", width=0.0005, position = position_dodge(width = 0.0009))+
    geom_smooth(method="lm")+
                theme_cowplot()+
  xlab("Climate Decomposition Index")+
  ylab("Soil Oxalate-Extractable Iron"
  )

ggarrange(basal_area_cdi, feox_cdi, ncol=2, nrow=1, labels= c("a","b"), legend=F)
## `geom_smooth()` using formula 'y ~ x'

ggsave("FigS2.pdf", path="figures/", width= 24, height= 12, units="cm")

summary(lm(plotBA~CDI,data=final_data )) #
## 
## Call:
## lm(formula = plotBA ~ CDI, data = final_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.7455 -0.7491 -0.1379  0.6662  2.1157 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.7240     0.4667   3.694 0.000596 ***
## CDI           0.5506    12.8459   0.043 0.966001    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9489 on 45 degrees of freedom
## Multiple R-squared:  4.082e-05,  Adjusted R-squared:  -0.02218 
## F-statistic: 0.001837 on 1 and 45 DF,  p-value: 0.966
summary(lm(FeOx~CDI,data=final_data )) #yes, more FeOx in warmer sites
## 
## Call:
## lm(formula = FeOx ~ CDI, data = final_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.5624 -2.8296  0.2597  1.6610 11.9572 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -1.093      1.821  -0.600 0.551306    
## CDI          184.539     50.119   3.682 0.000617 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.702 on 45 degrees of freedom
## Multiple R-squared:  0.2315, Adjusted R-squared:  0.2144 
## F-statistic: 13.56 on 1 and 45 DF,  p-value: 0.0006175
summary(lm(mg.C.per.g.soil_HF~FeOx,data=final_data %>% filter(Site=="LENO") )) #yes within this site Fe Ox is a big driver of concentrations
## 
## Call:
## lm(formula = mg.C.per.g.soil_HF ~ FeOx, data = final_data %>% 
##     filter(Site == "LENO"))
## 
## Residuals:
##       1       2       3       4       5       6       7 
## -4.0951  4.7057 -1.6980  2.2892 -1.6658  0.1938  0.2702 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   3.5803     2.5703   1.393  0.22239   
## FeOx          0.9075     0.2157   4.207  0.00843 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.16 on 5 degrees of freedom
## Multiple R-squared:  0.7798, Adjusted R-squared:  0.7357 
## F-statistic:  17.7 on 1 and 5 DF,  p-value: 0.008429
summary(lm(mg.C.per.g.soil_HF~FeOx,data=final_data %>% filter(Site=="DELA") ))  #yes within this site Fe Ox is a big driver of concentrations
## 
## Call:
## lm(formula = mg.C.per.g.soil_HF ~ FeOx, data = final_data %>% 
##     filter(Site == "DELA"))
## 
## Residuals:
##       1       2       3       4       5       6       7 
## -1.8847 -0.2968  3.2730  3.2409 -4.2195  1.1282 -1.2411 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   2.6878     4.0998   0.656   0.5410  
## FeOx          1.1080     0.4278   2.590   0.0488 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.016 on 5 degrees of freedom
## Multiple R-squared:  0.5729, Adjusted R-squared:  0.4875 
## F-statistic: 6.707 on 1 and 5 DF,  p-value: 0.04884
summary(lm(FeOx~E,data=final_data )) #no overall trend with FeOx and ECM dominance across sites
## 
## Call:
## lm(formula = FeOx ~ E, data = final_data)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -5.211 -2.730 -0.970  2.103 13.966 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    4.182      1.125   3.716 0.000557 ***
## E              2.301      1.933   1.190 0.240214    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.158 on 45 degrees of freedom
## Multiple R-squared:  0.03052,    Adjusted R-squared:  0.008974 
## F-statistic: 1.417 on 1 and 45 DF,  p-value: 0.2402
summary(lm(FeOx~E,data=final_data %>% filter(Site=="LENO") )) # yes within LENO positive relationship between soil FeOx and ECM dominance
## 
## Call:
## lm(formula = FeOx ~ E, data = final_data %>% filter(Site == "LENO"))
## 
## Residuals:
##       1       2       3       4       5       6       7 
##  6.4818 -0.3810 -0.9449 -2.4064 -1.4924  3.2959 -4.5530 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)    4.001      2.782   1.438   0.2099  
## E             11.881      4.206   2.825   0.0369 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.066 on 5 degrees of freedom
## Multiple R-squared:  0.6148, Adjusted R-squared:  0.5378 
## F-statistic: 7.981 on 1 and 5 DF,  p-value: 0.03689
summary(lm(FeOx~E,data=final_data %>% filter(Site=="DELA") )) # not at all within DELA 
## 
## Call:
## lm(formula = FeOx ~ E, data = final_data %>% filter(Site == "DELA"))
## 
## Residuals:
##       1       2       3       4       5       6       7 
## -5.2936 -0.5310 -0.5159 -0.7887  1.2474  2.2708  3.6111 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)    8.768      1.978   4.433  0.00681 **
## E              1.162      4.219   0.275  0.79401   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.129 on 5 degrees of freedom
## Multiple R-squared:  0.01495,    Adjusted R-squared:  -0.1821 
## F-statistic: 0.07586 on 1 and 5 DF,  p-value: 0.794
#Site means of MAOM C and N concentrations and proportions, ordered by climate decomposition index

MAOM_C_conc_cdi=ggplot(final_data, aes(x=CDI, y=mg.C.per.g.soil_HF))+
   geom_point()+
  stat_summary(fun= mean, fun.min=mean, fun.max=mean, geom="crossbar", width=0.0009, position="dodge")+
  stat_summary(fun.data = mean_se, geom = "errorbar", width=0.0005, position = position_dodge(width = 0.0009))+
  geom_smooth(method="lm")+
                theme_cowplot()+
  xlab("Climate Decomposition Index")+
  ylab("[MAOM C] (mg/g soil)"
  )

MAOM_C_prop_cdi=ggplot(final_data, aes(x=CDI, y=proportion.C_HF))+
   geom_point()+
  stat_summary(fun= mean, fun.min=mean, fun.max=mean, geom="crossbar", width=0.0009, position="dodge")+
  stat_summary(fun.data = mean_se, geom = "errorbar", width=0.0005, position = position_dodge(width = 0.0009))+
  theme_cowplot()+
  xlab("Climate Decomposition Index")+
  ylab(expression("proportion C"[MAOM])
  )

MAOM_N_conc_cdi=ggplot(final_data, aes(x=CDI, y=mg.N.per.g.soil_HF))+
   geom_point()+
  stat_summary(fun= mean, fun.min=mean, fun.max=mean, geom="crossbar", width=0.0009, position="dodge")+
  stat_summary(fun.data = mean_se, geom = "errorbar", width=0.0005, position = position_dodge(width = 0.0009))+
  theme_cowplot()+
  xlab("Climate Decomposition Index")+
  ylab("[MAOM N] (mg/g soil)"
  )

MAOM_N_prop_cdi=ggplot(final_data, aes(x=CDI, y=proportion.N_HF))+
   geom_point()+
  stat_summary(fun= mean, fun.min=mean, fun.max=mean, geom="crossbar", width=0.0009, position="dodge")+
  stat_summary(fun.data = mean_se, geom = "errorbar", width=0.0005, position = position_dodge(width = 0.0009))+
  geom_smooth(method="lm")+
                theme_cowplot()+
  xlab("Climate Decomposition Index")+
  ylab(expression("proportion N"[MAOM])
  )

ggarrange(MAOM_C_conc_cdi, MAOM_C_prop_cdi, MAOM_N_conc_cdi, MAOM_N_prop_cdi, ncol=2, nrow=2, labels= c("a","b", "c", "d"), legend=F)
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

ggsave("FigXXX.jpg", path="figures/", width= 24, height= 22, units="cm")
#Is the amount of oxalate-extractable iron in soil related to mycorrhizal type?

ggplot(final_data, aes(x=E, y=FeOx))+
  geom_point(color="blue")+
  geom_smooth(method="lm", color="blue")+
  facet_wrap(~Site)+
  theme_cowplot()
## `geom_smooth()` using formula 'y ~ x'

ggsave("ECM_vs_FeOx.jpg", path="figures/", width= 30, height= 24, units="cm")
## `geom_smooth()` using formula 'y ~ x'
ggplot(final_data, aes(x=CDI, y=FeOx))+
  geom_point(color="blue")+
  geom_smooth(method="lm", color="blue")+
  theme_cowplot()
## `geom_smooth()` using formula 'y ~ x'

ggsave("CDI_vs_FeOx.pdf", path="figures/", width= 30, height= 24, units="cm")
## `geom_smooth()` using formula 'y ~ x'
  C.N_FLF C.N_OLF C.N_HF
Predictors Estimates std. Error std. Beta standardized std. Error p df Estimates std. Error std. Beta standardized std. Error p df Estimates std. Error std. Beta standardized std. Error p df
(Intercept) 23.724 3.861 0.000 0.174 <0.001 41.000 22.582 5.534 -0.001 0.199 <0.001 41.000 23.195 2.405 -0.001 0.160 <0.001 41.000
E 7.898 2.508 0.413 0.131 0.003 41.000 12.146 3.133 0.494 0.127 <0.001 41.000 2.322 1.110 0.171 0.082 0.043 41.000
CDI 66.397 108.571 0.119 0.195 0.544 41.000 141.022 156.075 0.197 0.218 0.372 41.000 -299.874 67.826 -0.757 0.171 <0.001 41.000
FeOx -0.641 0.245 -0.441 0.168 0.012 41.000 -0.428 0.319 -0.229 0.171 0.187 41.000 -0.013 0.118 -0.013 0.114 0.912 41.000
Random Effects
σ2 25.801 39.696 4.909
τ00 3.948 Site 10.891 Site 2.607 Site
ICC 0.133 0.215 0.347
N 7 Site 7 Site 7 Site
Observations 47 47 47
Marginal R2 / Conditional R2 0.241 / 0.341 0.220 / 0.388 0.611 / 0.746
  C.N_HF C.N_HF C.N_HF C.N_HF C.N_HF C.N_HF C.N_HF
Predictors Estimates p Estimates p Estimates p Estimates p Estimates p Estimates p Estimates p
(Intercept) 20.25 0.001 16.55 0.002 13.96 0.003 10.98 <0.001 10.58 0.006 10.00 <0.001 9.77 <0.001
E -2.03 0.701 3.25 0.347 4.18 0.299 3.25 0.221 4.71 0.331 -0.46 0.540 -0.76 0.183
Observations 7 6 6 7 7 7 7
R2 / R2 adjusted 0.032 / -0.162 0.221 / 0.026 0.262 / 0.078 0.281 / 0.138 0.188 / 0.025 0.080 / -0.104 0.323 / 0.188

`